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Abstract 



We establish inversion formulas of the so called filtered back-projection type to 
recover a function supported in the ball in even dimensions from its spherical means 
over spheres centered on the boundary of the ball. We also find several formulas to 
recover initial data of the from (/, 0) (or (0, g)) for the free space wave equation in 
even dimensions from the trace of the solution on the boundary of the ball, provided 
the initial data has support in the ball. 
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1 Introduction and Statement of Results 



The problem of determining a function from a subset of its spherical means has a rich 
history in pure and applied mathematics. Our interest in the subject was provoked by the 
new medical imaging technologies called thermoacoustic and photoacoustic tomography. The 
idea behind these [TU1 [16] is to illuminate an object by a short burst of radiofrequency or 
optical energy which causes rapid (though small in magnitude) thermal expansion which 
generates an acoustic wave. The acoustic wave can be measured on the periphery or in 
the exterior of the object. The inverse problem we consider is to find the distribution of 
the absorbed energy throughout the body. This is of interest, since the amount of energy 
absorbed at different points may be diagnostic of disease or indicative of uptake of probes 
tagged to metabolic processes or gene expression For a more thorough discussion of 
the modelling and biomedical applications, the reader is referred to the recent survey |17j . 
If the illuminating energy is impulsive in time, the propagation may be modelled as an 
initial value problem for the wave equation. The problem of recovering the initial data of a 
solution of the wave equation from the value of the solution on the boundary of a domain is of 
mathematical interest in every dimension, but for the application to thermo-/photoacoustic 
tomography it would appear that the three dimensional case is the only one of interest, since 
sound propagation is not confined to a lower dimensional submanifold. However, there exist 
methods of measuring the generated wave field which do not rely on point measurements 
of the sort that would be generated by an (idealized) acoustic transducer. In particular, 
integrating line detectors, which have been studied in [3J [H] , in effect compute the integral 
of the acoustic wave field along a specified line. In this paper, we work under the assumption 
that the speed of sound, c, is constant throughout the body, and since the x-ray transform 
in a given direction of a solution of the three dimensional wave equation is a solution of 
the two dimensional wave equation, the problem is transformed. If a circular array of line 
detectors is rotated around an axis orthogonal to the direction of the line detectors [?1 IT3] , 
then for each fixed rotation angle the measurement provides the trace of the solution of the 
two dimensional wave equation on the circle corresponding to the array. The initial data of 
this two dimensional problem is the x-ray transform of the three dimensional initial data. If 
the inital data can be recovered in the disk bounded by the detector array and assuming that 
the projection of the object to be imaged lies in this disk, then the problem of recovering 
the three dimensional initial data is reduced to the inversion of the x-ray transform in each 
plane orthogonal to the axis of rotation. One such two dimensional problem is illustrated in 
Figure dJ 

To our knowledge, the first work to tackle the problem of recovering a function from its 
circular means with centers on a circle was [IB], whose author was interested in ultrasound 
reflectivity tomography. He found an inversion method based on harmonic decomposition 
and for each harmonic, the inversion of a Hankel transform. This method has been the basis 
for most subsequent work on exact inversion of circular means. The inversion of the Hankel 
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Figure 1: Principle of thermoacoustic tomography with integrating line detectors. 

A cylindrical array of line detectors records the acoustic field and is rotated around the axis 
e±. For fixed rotation angle the array outputs the x-ray transform (projection along straight 
lines) of the solution of the wave equation restricted to lines passing through the boundary 
S of the disk. The initial condition is given by the x-ray transform of the initially induced 
pressure restricted to lines orthogonal to the base of the cylinder. 

transform involves a quotient of a Hankel transform of a harmonic component of the data 
and a Bessel function. That this quotient be well-defined turns out to be a condition on the 
range of the circular mean transform [2] . See also [1] for range results on the spherical mean 
transform on functions supported in a ball in all dimensions, and \6\ for range results for the 
wave trace map for functions supported in the ball in odd dimensions. 

In the work of the first and third authors with Sarah Patch plj, several formulas were 
found to recover a smooth function / with support in the closure B of the open ball B C R n 
from the trace of the solution of the wave equation on the product dB x [0, diam(i?)] provided 
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that the space dimension is odd. Specifically, if u is the solution of the initial value problem 



u tt - Au = 0, in R n x [0, oo) (1) 

u(.,t = 0) = /(.), u t (.,t = 0)=0, (2) 

where / is smooth and has support in the B, then several formulas were found to recover / 
from u(p, t) for p E S := dB and t G R + . 



The first and third authors tried, at that time, to extend the method to even dimensions, 
but did not see a way. Recently, the second author tried numerical experiments using a 
two dimensional analog of one of the inversion formulas and found that it gave excellent 
reconstructions. This prompted our re-examination of the problem. Among the results of 
this paper is a proof of the validity of this formula. 



To describe our results, we introduce some notation. The spherical mean transform Ai 
is defined by 

(Mf)(x,r) = -^r I f(x + re)dS(e) (3) 
l D I J s n - x 

for / G C°°(R n ) and (x,r) G R™ x [0, oo). In this expression, |S' n_1 | denotes the area of the 
unit sphere S*™ -1 in R n and dS(6) denotes area measure on the sphere. In general, we write 
the area measure on a sphere of any radius as dS, except when n = 2 when we write ds. We 
will denote the (partial) derivative of a function q with respect to a variable r by d r q, except 
in a few formulas where the subscript notation q r is used. At several points we use D r to 
denote the operator 

(D r u)(r) := MM 

acting on smooth (even) functions u with compact support. Moreover, r will be used to 
denote the operator that multiplies a function u(r) by r. 



Our first set of results is a pair of inversion formulas for the spherical mean transform in 
even dimensions. We state and prove these first in dimension two; that is, for the circular 
mean transform. 

Theorem 1. Let D C R 2 be the disk of radius R centered at the origin, let S := 3D denote 
the boundary circle, and let f G C°°(R 2 ) with supp / C D. Then, for x G D, 

f(x) = — -A x / r(Mf)(p,r)\og\r 2 ~\x-p\ 2 \ drds(p) (4) 

and 

f( x ) = 7T-B- / (d r rd r Mf)(p,r)log\r 2 - \x - p\ 2 \ drds(p). (5) 
2tt-Ko Js Jo 
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In Theorem [H d r rd r Aif denotes the composition of d r , r, d r and M. applied to /. 
The same convention will be used throughout the article to denote the composition of any 
operators. 



While A4 f has a natural extension to the negative reals as an even function, we instead 
take the odd extension in the second variable. Then formula (j5]) has the following corollary: 

Corollary 1. With the same hypotheses as in Theorem^ and M. f extended as an odd 
function in the second variable r, f can be recovered for x G D by 

m i /f!«„ s(p)i (6) 

2ttR J s J_ 2Ro \x-p\-r 

and 

ff \ 1 f\ I f 2R ° ( d rMf)(p,r) j 

f( x ) = 7T-5- / \x-p\ / — j ; drds{p), (7) 



2tiR J s J_ 2Ro \x-p\-r 

where the inner integrals are taken in the principal value sense. 

These forms are very close to the standard inversion formula for the Radon transform in 
the plane [32J Eq. (2.5)]. 

In higher even dimensions we prove a similar pair of results. 

Theorem 2. Let B C R n , n > 2 even, be the ball of radius Rq centered at the origin, let 
S := dB be the boundary of the ball, set 

Cn = (-l)("- 2 )/ 2 2((n - 2)/2)!7r"/ 2 = (-l)( n - 2 '/ 2 [((n - 2)/2)!] 2 |5 n - 1 |, 

and let f G C°°(R n ) have support in B. Then, for x G B, 
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2R 



f(x) — — — A x / / log|r 2 - \x-p\ z \ (rD?- 2 r n - 2 Mf)(p,r)drdS(p), 
c n Ro Js Jo 
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JS JO 



2R 

f{x) = -=— / / \og\r 2 -\x-p\ 2 \{rD r ;r 1 r n - 1 d r Mf){p,r)drdS{p). (9) 



Recently, Kunyansky [TT] has also established inversion formulas of the filtered back- 
projection type for the spherical mean transform. His method and results appear to be very 
different than ours. 

For some results, it will be more convenient to use the wave equation (CQ) with initial 
condition 

M (.,t = 0) = 0, ut{.,t = 0) = /(.), (10) 
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It is obvious that the solution of §B) with initial values ([2} is the time derivative of the 
solution of (pQ) with initial values (fTUj) . We denote by V the operator which takes smooth 
initial data with support in B to the solution of ([1]), (flOj) restricted to S x [0, oo) and by W 
the operator taking / to the solution of ([T]), restricted to 5 X [0, oo). These operators are 
simply related by W = d t V . An explicit representation for V comes from the well-known 
formula [I] 



giving the solution of the initial value problem ([T]), (fTUj) . in dimension n > 2. We denote by 
V* and W* = —V*d t the formal L 2 adjoints of V and W mapping from smooth functions 
u G C°°(S x [0, oo)) with sufficient decay in the second variable. An explicit expression for 
V* u will be given in Section [3j 

We have two types of inversion results for the wave equation. The first type is based on 
the inversion results for the spherical mean transform, since the spherical mean transform 
itself can be recovered from the solution of the wave equation by solving an Abel type 
equation. In dimension two, this approach yields the following result. 

Theorem 3. Let D C R 2 be the open disc with radius Rq and let S := dD denote the 
boundary circle. Then there exists a kernel function K : [0,2Rq] 2 — > R such that for any 
f G C°°(R 2 ) with support in D and any x G D 



An analytic expression for K will be given in Section 

Theorem [3] provides inversion formulas of the filtered back-projection type for recon- 
struction of / from (W/)(p, t) = (d t V f)(p,t) using only data with t G [0,2_R ], despite the 
unbounded support of W / and V f in t. 

The second type of inversion results holds in all even dimensions and takes the following 
form. 

Theorem 4. Let f be smooth and supported in closure of the ball B of radius Rq in R 2m , 
and let V f and W / be as above. Then for x G B 




(11) 




(12) 



f(x) = -^(V*td 2 Vf) (x) 

/(*)= j-(W*tW/)(x): 



^(V*d t td t Vf) (x). 



(13) 



(14) 
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We will prove ({TBI in dimension n = 2m = 2 directly. The higher dimensional case of 
(fl3|) . and (jbil) in all dimensions, are consequences of the following trace identities, relating 
the L 2 inner product of the initial data to the weighted L 2 inner product of the traces of the 
solutions of the wave equation. 

Theorem 5. Let f,g be smooth and supported in the ball B of radius Rq, in R 2m with 
m > 1, let S := OB, and let u (resp. v) be the solution of the initial value problem (T7J), / TiOj) 
with initial value f (resp. g). Then 

[ f(x)g(x)dx = -— I [ tu u (p,t)v(p,t)dtdS(p), (15) 
Jb tt o Js Jo 

/ f(x)g(x)dx= -5- / / tu t (p,t)v t (p,t)dtdS(p). (16) 
Jb -^0 Js Jo 



In the proof of this theorem, (Tl5l) for n = 2 follows from ( Tl3l) for n = 2, while ( T151) in 
higher even dimensions is derived from the n = 2 case; (TlBl) is a consequence of (fl5l) in all 
dimensions. We remark that these identities were already proved in [5] for odd dimensions, 
and so they hold for all dimensions. 

Section 2 is devoted to the proof of the inversion formulas for the spherical mean trans- 
form, that is, Theorems 1, 2, and Corollary 1. Section 3 treats the wave equation and 
contains the proofs of Theorems 3, 4, and 5. This is followed by a section reporting on the 
implementation of the various reconstruction formulas of the preceding sections and results 
of numerical tests, in dimension two. 



2 Spherical Means 



In this section we prove the Theorems related to the inversion from spherical means and 
Corollary [lj We begin by establishing an elementary integral identity, which is the key to 
the results in this paper. 

Proposition 2.1. Let D C R 2 be the disk of radius Rq, and let S = dD be the boundary 
circle. Then for x, y G D with x ^ y, 

log I \x — p\ 2 — \y — p\ 2 \ ds(p) = 2TcR \og \x — y \ + 2TcR \ogR . (17) 



Proof. Let x ^ y both lie in D and let I denote the integral on the left on (fl71) . Expanding 
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the argument of the logarithm as 



\x — p\ — \y — p\ 



2R \x - y\ 



x + y p \ x — y 



2R R J \x-y 



setting e := ifz^p an d writing p = R 9 for 6 G S* 1 , we have 



/ = 2tt.Ro log (2R \x -y\) + Ro log |e • 9 - a\ d8, 

Js 1 



where 



x + y 
a = • e 

2Rn 



\x 


2 




y 


2 


2R 


x - 


-y\ 



We note that \a\ < 1. 

Using the parameterization 6 = cos(0)e + sin(0)e _L , the integral term on the right of ( fl8l) 
has the form 

p2tt 

Ro / log I cos <fi — a\ dtp. 
Jo 

Writing a = cos a and using the sum to product trigonometric identity cos0 — cos a = 
—2 sin ((0 + a)/2) sin ((0 — a)/2), this is equal to 

i?o / (log 2 + log | sin ((0 + a)/2) | + log | sin ((0 - a)/2) |) 
Jo 

By periodicity, and two linear changes of variable, this reduces to 

p2n PIT 

R / (log 2 + 2 log | sin (0/2) |) d(j) = 2R n log 2 + AR j log sin u du, 
Jo Jo 

which is independent of a, and hence of x and y. The latter integral in can be found in 
tables, and is equal to — Roii\og2 ) so the sum is —2tcR log 2. Substituting in ([IB]) gives the 
desired result. □ 



Proposition 12.11 is already enough to establish Theorem [TJ 

Proof of Theorem [TJ. Let / G C°°(R 2 ) be supported in D and let p be any point in 
S = dB. Using the definition of M. f and Fubini's theorem, we have that 



2R 



(r M. /) (p, r)q(r) dr 



2tt 



f(p + z)q(\z\)dz, 



(19) 



R 2 
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for any measurable function q provided that the product of functions on the right is absolutely 
integrable. Applying this with q(r) = log |r 2 — \x — p\ 2 \ and making the change of variables 
y = p + z gives 

(r j\4)(f)(p, r) log |r 2 — \x — p\ 2 \ dr ds(p) 

= ^~ / f(y) l °g\\y-p\ 2 - \ x -P\ 2 \ dyds{p). 
^ Js Jri 

Fubini's theorem again justifies the change of order of integration in the iterated integral on 
the right hand side, and so 

1 f /" t 1 1 _ . _.|2 I 121 j.. 2nR 



7T / /(y) / !og - |ac - p| tZs(p) dy = — — / /(y)(log \x - y\ + \ogR ) dy 

^ J Js 27r Jr. 2 

upon application of ffTTj) . Recalling that for any constant c, l/(27r) log |x — y| + c is a 
fundamental solution of the Laplacian in R 2 , we have 
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2R 



f( x ) = 7T^r A x / / (rMf)(p,r)\og\r - \x - p\ \drds(p), 
2irK Js Jo 

which proves (j3J). 

The second formula, (131) , has a similar proof. In this case, we use that the spherical 
means satisfy the Euler-Poisson-Darboux equation [1] 

(d 2 r Mf)(x,r) + - (d r Mf)(x,r) = (AM f) (x,r) = (M Af)(x,r). 
r 

The left hand side of the Darboux equation may be written as (l/r)(d r rd r Ai f)(x,r), so the 
expression on the right of (jSJ) may be rewritten as 

[r M. Af)(j>, r) log |r 2 — \x — p\ 2 \ drds(p). (20) 



2ttR 




o Js Jo 



Again applying (I19p . now with the function q(r) = r log |r 2 — \x — p\ 2 \ and Af instead of /, 
interchanging the order of integration and using (fTTj) shows that the expression (120]) is equal 
to 

7T / \f(y)( lo s\ x -y\ + lo sRo)dy = f(x), 

2vr J R2 

since no boundary terms arise in view of the support hypothesis on /. □ 
Proof of Corollary [TJ Let x 6 D, and let 

f 2Ro 

U(p,x) := / (d r rd r M. f) (p, r) log |r 2 — |x — p| 2 1 o?r 



denote the inner integral in (JSj). Taking the support of / into account, writing the logarithm 

as 

log |r 2 — \x — p\ 2 \ = log \r — \x — p\ | + log |r + \x — p\\ , 
and integrating by parts leads to 

um = -p.v. f^M»- r(ra r Mf)(pr) dr 

Jo r-\x-p\ J r + \x-p\ 

Here we have used that the distributional derivative of log |r| is P.V. 1 as well as an ordinary 
integration by parts. Therefore §5$) implies 

(21) 

where the inner integral of the first term on the right is taken in the principal value sense. 
The odd extension of A4 f, M. f(p, — r) := — M. f(p, r), is smooth on R since M. f vanishes 
to infinite order at r = by the support hypothesis on / and (rd r A4 f)(p,r) is an odd 
function in r. Substituting r = —r in the second integral in (j2ip gives 

f(x) = J±. f p(^^/)far) drd + ^Lf f° (ra r Mf)M 
2-kRqJsJq r-\x-p\ 2ttR J s J_ 2Ro r - \x - p\ 

and hence 

2ttR J s J_ 2Ro \x-p\-r 
This is dSD- To prove ([7]), it suffices to write 

\x — p\ 



-1 + 



\x — p\ — r |x — p\ — r 



in ([6]) and to note that /^ (<9 r A4 f)(p, r) dr = 0, by the support hypothesis on /. □ 



2.1 Proof of Theorem [2 



We have found several proofs of Theorem [2J the extension of Theorem [T] to higher even 
dimensions. The one we present is based on reduction of the higher dimensional problem to 
the two dimensional case already established. Another, which is not presented in this article, 
is based on an extension of (TlT|) to higher dimensions. 
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We first observe that by a dilation, we may reduce the problem to the case when / is 
supported in the unit ball. Tracing through the formulas (|HJ) and (Q it is routine to verify 
that scaling from the unit ball to the ball of radius Rq introduces a factor of Rq. To simplify 
notation, we shall now suppose that / is supported in the unit ball B. Let Q and iV denote 
the operators 

(Qf)(x) = A x [ [(rD?- 2 r n - 2 Mf)(p,r)\og\r 2 -\x-p\ 2 \drdS(p), (22) 
Js Jo 

(Nf){x) = [ [ (rD^ 1 r n ~ 1 d r Mf)(p,r)\og\r 2 ~\x~p\ 2 \drdS(p), (23) 
Js Jo 

that map / G C°°(R n ) supported in B to constant multiples of the the right hand sides of 
dHJ) and 0. Moreover (/, g) denotes the L 2 product of two functions supported in B. To 
establish Qf = c n f and Nf = (c n /2)f we will use the following auxiliary results. 

Proposition 2.2. Let f , g be smooth and supported in B. Then 

[ (Qf)(x)g(x)dx=(Qf,g) = 2(f,Ng) = 2 [ f(x)(Ng)(x)dx. (24) 



Proof. Let F = M. f and G = AA. g. Using the self-adjointness of A, applying Fubini's 
theorem and an n-dimensional analogue of f|T9|) . we obtain 

{Qf, g) = j B yj\rD?- 2 r n - 2 F)(p,r)log\r 2 -\x-p\ 2 \ drdS(p)^J (A x g)(x) dx 

= \S n ~ x \ f f o (f o (rD?- 2 r n - 2 F)(p,r)\og\r 2 -f 2 \(MA x g)(p,f)f n -Uf^J drdS{p) 

= IS™" 1 ] J jf Qf (rD?- 2 r n - 2 F)(p,r)log\r 2 - f 2 \ dr^j (M A x g)(p, f)f n_1 dfdS(p) 

= \S n - l \ f J q (J q (rD?- 2 r n - 2 F)(p,r)\og\r 2 - f 2 \ dr^j d f f n ~ l d f G{p,f) df dS{p). 

(25) 

To justify the last equation it is used that G satisfies the Euler-Poisson-Darboux equation and 
the identity f n ~ 1 (d 2 + It ^df) = df{f n ~ l df). Applying the identities (D™~ 2 )*r log \r 2 — f 2 \ = 
(— l) n ~ 2 rD™~ 2 log |r 2 — f 2 | = rDf~ 2 \og\r 2 — f 2 \ in two stages to the last expression, this 
becomes 

IS" -1 ! J J (J r n - x F(p,r)D^- 2 \og\r 2 -r 2 \dr\ (drf n - x drG{p,f)) dfdS{p) 
= J J 2 (Jj{y)D^ 2 \og\\y-p\ 2 -f 2 \dy^j (d P f n -%G(p,f)) dfdS{p) 
= J(J 3 f Q lo S I \V ~ P\ 2 ~ f2 \ m) n -%r n -%G)(p, f) drdS(jfi) f(y) dy 
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after applying Fubini's theorem. This is finally seen to be equal to (/, 2Ng) since (D™ )*<9 f = 
2f(-l) n - 2 D?-\ □ 



We now look at the spherical means of products 

f(x)=p k a(pM6), (26) 

where x = p9 with p > 0, 9 G S n ~ l , $ is a spherical harmonic of degree k, and a : R — > R 
is an even smooth function supported in [—1,1]. Let F := M. f be extended to an even 
function in the second component and let v — n + 2k. Then F satisfies the initial value 
problem (IVP) for the Euler-Poisson-Darboux equation 

77,-1 \ 

dlF + d r F (x, r) = A x F(x, r), (x, r) G R" x R (27) 

r J 

F(x, 0) = a(p)p k §(9), d r F(x, 0) = 0, x G R n , (28) 

and, conversely, any solution of ( 1271) . ( |28l) . is the spherical mean of the initial values. The 
unique solution of ( 1271) . ( 1281) has the form F(x,r) = p k A(p,r)§(9) where A(p,r) is the 
solution of the IVP 

(L n A)(p,r) = (d 2 p A+^d p A^ (p,r), (p,r) G R 2 , (29) 
A(p,0) = a(p), d p A(p,0) = 0, pGR. (30) 
Here (L n A)(p,r) := (d 2 A + ^d r A)(p,r). 

We recall that the operator D r satisfies L n D r = D r L n ^2 and for any p G N 



<9 2 + ^<9 r (r^) = r» d 2 r + —^d r w, 

r J \ r J 

that is L 2 _^r M = r At L M+2 . So 

(L 2 _ M+2CTJ D?r^)(r) = p'L^r"™) (r) = (D^L^ 2 w)(r). (31) 

If we set p = n — 2 and cr = (n — 2)/2 in (I3TI) . then p + 2 = n and 2 — p + 2cr = 2. Therefore 

(L 2 D r (n - 2 )/ 2 r n ~ 2 u;)(r) = (Z^ n - 2)/2 r n - 2 L,p/;)(r). (32) 



Now we set 



B(p,r) := i{n _ 2)/2 y ( D ^ 2)/2rn ' 2A )(P^)- (33) 
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Since A(p, r) is even in r and D r corresponds to differentiation with respect to r 2 , H(p, r) is 
even in r. Moreover, by (j30~l) . H(p,0) = ^ n _\y 2 y_ A(p, 0) (D r n 2 ^ 2 r n ~ 2 ) = a(p), and therefore 
from ([29]) and ([32]) it follows that H is the solution of the IVP 

d 2 r H + ~d r Hj (p, r) = (d 2 p H + ^d P H) (P, 0, (P, e R 2 , (34) 
F(p, 0) = a(p), d r H(p, 0) = 0, p G R, (35) 

Proposition 2.3. Lei Aj(p,r) ; i — 1,2 solve ( fffPj) n = 2, subject to initial conditions 
Aj(p, 0) = ctj(p), <9 r Aj(p, 0) = 0, where ctj are smooth even functions with support in [— 1, 1] 
and v >2 is even. Then 

1 ,2 ,2 

p" 1 «i(p)a 2 (p) dp = — / / rAi(l,r)dflog\r 2 — r 2 \r(dfA 2 )(l,r)drdr. (36) 

Jo io 



Proof. Let /c = (z/ — 2)/2 and let $(#) be a nontrivial real circular harmonic of degree k. 
Then Fi(x, r) : = Ai{p,r)p k <S>(6) satisfies (ETJ), (EH]) for n = 2, and so is the circular mean of 
its initial value, fi(x) = ai(p)p k <S?(6). By (TjJ, fi = ^Qfi, and using ([25jl gives 

(/i,/ 2 > = ^-Wi,/ 2 > 

= / / / rFi(p,r) log |r 2 — f 2 |(<9 f f<9 f F2)(p, f) drdf ds(p). 
Js Jo Jo 

Taking account the form of and that p = 1 on S, this may be rewritten as 

(f 1 ,f 2 )=[$ 2 (p)ds(p)[ [ rA!(l,r)log|r 2 - f 2 |(<9 f f<9 f A 2 )(l,f) drdf. (37) 
Appealing to the form of fi = F(x, 0), 

A ■ h a-,)dp [ $ 2 (p)ds(p) 

(38) 



(fi,h)= I p(p k a 1 )(p k a 2 )dp j^ 2 (p)ds(p) 
p u - 1 a 1 (p)a 2 (p)dp / $ 2 (p)ds(p). 



Js 

Since L $ 2 (p) ds(p) ^ 0, a comparison of (l37j) and (1381) and an integration by parts on the 
right side of (1371) establishes (1361) which completes the proof. □ 

Proof of Theorem [21 Let {$•,} be an orthonormal basis for the spherical harmonics on 
S n ~ l , and consider fi, i — 1, 2 of the form (1261) with a = q.{ and $ = $^ of possibly different 
degrees. Let Fi be the even extensions of M. fi as above. Then by orthogonality, (fx, f 2 ) = 
unless ji = j 2 , in which case 

(A,/ 2 >= ^ pf , - 1 a 1 (p)a 2 (p)dp, (39) 
./o 
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with v = n + 2k, where k is the degree of Evaluating (Qfi, / 2 ) by (1231) . and using that 
Fi = p ki Ai(p,r)Qj t , we see that it is also zero unless j\ = j%. In this case we have 

(Q/1,/2) = IS"" 1 ! f [\rD^ 2 r n - 2 A 1 )(l,r)log\r 2 -f 2 \(d^-%A 2 )(l,f)drdf 
Jo Jo 

= \S n - 1 \ [ [ {D^r n - 2 A 1 ){l,r){D;)^(rlog\r 2 -f 2 \)(drf^- 1 d f A 2 )(l,f)drdf 
Jo Jo 

2 /" 2 ra-2 rt-2 

n— 1| / I „Cn 2 „n-2 a \/i „\ n 2 1 l„2 ^2|/o=n— 1; 



= 1-5 1/ / r(D r 2 r n - 2 Ai)(l,r)D f a log |r 2 - f i |(9 f f n - 1 9 f A 2 )(l, f) drdf, 

(40) 

since (Z)£™ 2 - )//2 )*r log |r 2 — f 2 | = (— \Y n ~ 2 ^ 2 rD'^ 2 log |r 2 — f 2 | — rD^ n 2 ^ 2 log |r 2 — f 2 |. Ap- 
plying the adjoint (distributional derivative) again in (j40l) . 



(Q/1,/2) = |S n-1 | / / r(D r 2 r n ~ 2 A 1 )(l,r)log|r 2 -f 2 |(,D r - 2 )*(9 f f"- 1 a r -A 2 )(l, f) dr df 
</o Jo 

n 2 n-2 „_ 9 n-2 

r(ZV r n - 2 A 1 )(l,r)log|r 2 -f 2 |(-l) — (<9 f D f 2 f^d^Xl, f) dr df 

r 2 r 2 n-2 n-2 

= I^K-l)"/ 2 / / r(ZV r n - 2 A 1 )(l,r)9 f log|r 2 -f 2 |( J D f 2 f n - 1 9 f A 2 )(l, f) dr dr. 
Jo Jo 

We now use the following identity, which is readily proved by induction, 

Di n - 2)/2 f n ~%q = fd f Dt 2)/2 f n - 2 q 
taking q = A 2 , and observe that defining Hi by (1331) with A = Ai, it then holds that 

(Q/1,/2) = |5 n - 1 |(-l) n/2 [((n-2)/2)!] 2 T l\H 1 {l y r)d f \og\r 2 -f 2 \f{d f H 2 ){l,f)drdf. 

Jo Jo 

By ( 134|) and (1331) the if, satisfy the hypotheses of Proposition 12.31 with initial data a,, so by 
fl36|) the expression on the right is equal to 

(-l)(«-W2| 5 n-l| [((n _ 2)/2)!] 2 f p ^ ai{p)a2{p)dp . 



Thus we have proved that for /j of the form above 

<Q/i,/ a > = (-l) ( "- 2)/2 |S- 1 |[((n-2)/2)!] 2 (/ 1 ,/ 2 ). (41) 

We note that the constant on the right is c n of Theorem [2j By linearity and orthogonality 
of spherical harmonics, this still holds when either f\ or / 2 is replaced by a finite linear 
combination of such functions. The set of finite linear combinations of functions of form 
(123]) is dense in L 2 , and so we have Qf = c n f in L 2 when / is a finite linear combination 
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of functions of the form ff26]) . Now, let g be smooth with with support in the unit ball. 
Applying Proposition 12.21 it follows that 

(f,Ng) = (l/2){Qf,g) = (c n /2){f,g). (42) 

for all / as above. Since (l4~2l) holds for a dense subset of functions / in L 2 (B), it implies 
that Ng = (c n /2)g almost everywhere in B. However, Ng is easily seen to be a continuous 
function, and so Ng = (c n /2)g holds pointwise in B, which is But if N is a multiple of 
the identity, then so Q, and the proof is complete. □ 



3 The Wave Equation 



We begin the analysis of recovery of initial data from the trace of the solution of the wave 
equation on the lateral boundary of the cylinder. As mentioned in the introduction, we 
have two types of inversion results. The first, Theorem [3j is really a corollary of one of the 
inversion formulas for circular means from the previous section. 

Proof of Theorem [31 Let u(x,t) to be the solution of the IVP (00), (j2J) in dimension two. 
Then by (TITJ), 

for p e S. We can recover the circular means from u by the standard method of inverting an 
Abel type equation. The details are not hard and may be found, for example, in [12]. The 
result is 

(Mf)(p,t)=* [j^L-dt. (43) 
n Jo Vr 2 - t 2 

Inserting (T43"]) into the inversion formula (j4]) for M. and applying Fubini's theorem, gives, for 

x e D, 



l 



2-Ro 



f(x) = A / / (r M. f)(p, r) log r — \x — p\ \drds(p). 

2nR J S J 

i . r r 2Ro r u( P ,t) , , 2 l2l , , 7 . , 

A / / r - log r — \x — p\ at ar asypj 



Ron 2 JsJo Jo Vr 2 - 1 2 



— -A / / u(p,t) / , log \r 2 — \x — p| 2 | dr dtds(p) 

JT 2 JsJo Jt V^t^ ' ' 



Ro^ 2 
1 

R 7V< 



■2Ro 

All u(p,t)K(t,\x-p\)dtds(p). 




S JO 
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Since u(p,t) = (W f)(p,t), this is (p2D, with 

K(t,f) := / Wlr 2 -r 2 \dr (44) 

The integral in (|44|) can be evaluated exactly. For the sake of completeness, we give the 
analytic expression. If we substitute r = ^Jt 2 + £ 2 in (|4"4"|) . then <ir = (£/r)<i£ and thus 



^(*,r)= / log |e 2 + (t 2 -f 2 ) | di 

Jo 

= \J 4_Rq — t 2 (-2 + log |4i2jj - f 2 |) + r(t, f), 



where 



, yr 1 — t 1 log v , = — t < r, 



2A/t 2 — r 2 arctan \ 1 4 t > r 



t 2_ f 2 

□ 



For the second type of inversion formula, we start by deriving a representation of the 
formal adjoint V*, for n = 2. For any continuous function G(p,t) on S x [0, oo) that has a 
small amount of decay as t — > oo, by Fubini's theorem, we have 



p poo 

(Vf,G)= (Vf)(p,t)G(p,t)dtds(p) 
Js Jo 



where 



1 / / G '° 

Js Jo 

1 f f°° ( f* f f(p + ru) 




G(p,t)[ / f(p + ru)ds(u)dr)dtds(p) 

27r Js Jo \Jo Vt 2 - r 2 J si J 

rdrdS{uj) I G(p,t) dt ds(p) 



27r Js Jo \Jo Js^ v 7 ^ 3 

i r r ( f f(v) 



[ f{v) X ({\y -p\< t}) dy ) G(p, t) dt ds(p) 

Jr 2 \ t 2 - \y - v\ 2 / 



2?r Js Jo Wr 2 \/t 2 -\y -p 



2tt Ir2^ V ^ {^Isl\y. 



G{PJ) dtds(p) I /(//),//, 



-Pi - 12/ -P 



G) (y) := / /" |2 eft <fe(p) • (45) 

The integral in (|45|) will be absolutely convergent for continuous G provided that G has a 
small amount of decay as t — > oo, for example if G(p,t) = 0(l/t a ), as t — > oo, for some 
a > 0. 
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Next, we note a differentiation formula for the fractional integral appearing in (TlTT) . 
Proposition 3.1. Let h be differentiable on [0,oo). Then, fort > 0, 

Jo t Jo 



Proof. Making the change of variable r = t£ in the integral on the left we have to evaluate 

d t f t=th{tOdt. 

Jo V 1- ? 



Here differentiation under the integral yields J Q (^C^'(^C) + h(t^)) which is equal to 

the expression on the right side after changing back to integration with respect to r = t£. □ 

Proof of (11311 in Theorem [4] for n = 2. We compute (V* td 2 V f)(x) for smooth / 
supported in B and x G B. The function td 2 V f has decay of order 1/t 2 as t — > oo, and so 
lies in the domain of V*. Using the definitions of V and V* , and relation ()46p . 



{V*td 2 t Vf){x) 



1 f f°° d 2 ( [ r (Mf)(p,r) dr \ tdtds(p) 



2nJsJ\x-p\ \Jo Vt 2 - r 2 J ^Jt 2 - \x -p 12 



1 f r d (l [ t r ( 9 rrMf)(p,r) dr \ tdtds(p) 



2lT JsJ\x-p\ \t Jo Vt 2 - r 2 J ^t 2 -\x-p 



Carrying out the differentiation in t using the chain rule, using again ( l46|) . and combining 
terms, the last integral can be rewritten as 

1 f f°° ( f r r (d r rd r r M f) (p, r) - r (d r r M f) ( p, r) ^ . ^ ^ ( ^ ; 



2tt 




t-^/t 2 — \x — p\ 2 \/t 2 — r 2 



'S^li-pl \Jo 

Using the identity 

d r rd r rh — d r rh = d r r(d r rh — h) = d r rrd r h = d r r 2 d r h 
and applying Fubini's theorem (V* td 2 V f)(x) is in turn is equal to 

i- f r r (d r r 2 d r Mf)(p,r) (r — ==JL=—== ) drds(p). 

l7[ JS J0 \Jmax(\x-p\,r) t y/ V — \X — p\ z \Jt l — r l 

The inner integral evaluates to 

TTT~ 1 log 

2r\x — p\ 
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r + 


x — p\ 


r — 


x — p\\ 



giving 



(V*tdfVf)(x) = - 



S \JQ 



(d r r 2 d r M. f) (p, r) log 



r + 


x — p 


r — 


x — p\\ 



dr 



ds{p) 
\x — p| 



(47) 



Treating the inner integral in principal value sense, and integrating by parts, it is equal to 
the limit as e — > of boundary terms 



{r 2 d r M. /)(p, r) log 



r + \x — p\ 
\x — p\ — r 



\x-p\-e 



+ 



(r 2 d r M f)(p, r) log 



plus the term 



(rd r M. f) (p, r) rd r log 



r + \x — p| 
r — \x — p\ 



dr. 



\x-y\+e 



R+\[|x-p|-e,|a;-p|+e] 



r + 


x-y\ 


r — 


x — p 



Using that .M / is smooth, flat at r = 0, and of bounded support in (0, oo), the limit of the 
boundary terms is zero. Using the identity 



rd r log 



r + 


x — p 


r — 


x — p 



|x — p\d r log |r 2 — |x — p| 



followed by another integration by parts, yields the sum of another pair of boundary terms 
and 



I t = —\x — p\ 



(d r rd r M. f) (p, r) log |r 2 — \x — p| 2 | dr. 



' R+\[\x-p\-e,\x-p\+e] 

The boundary terms again evaluate to zero as e — * while the integral I e converges to 

poo 

— \x—p\ I (d r rd r A4 f) (p, r) log I r 2 — \x — p| 2 1 dr. 



Inserting this into ff47|) and taking into account the support of M. f, gives 

-2R 



(V*tdfVf)(x) 



1 

4:71 




(d r rd r M. f) (p, r) log r — |x — p| rfr<is(p). 



5 jo 



In view of ([5]) of Theorem [H ffl3l) in Theorem H] is proved, for n = 2. 



□ 



Proof of Theorem [5l Formula ffl5l) . for = 2, an easy corollary of the result just es- 
tablished. Indeed, for /, g smooth with compact support in the closed disk of radius R , 
then 

(/. *> = -4- ( v * td " Vf,g) = - J- <^ 2 P /, P ^> , (48) 

Kq Ho 
which is fTTS]) for = 2, due to the definition of the operator V. 
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In ( TT5|) . the left hand side is symmetric in / and g, while the right side is not. Thus there 
is a companion identity, reversing the roles of u and v on the right. Taking the difference 
gives the equation 



= / / t(u tt v -uv tt )dtds(p). 
Js Jo 

Integrating by parts (the boundary terms vanish) yields 



= / / (u t v — uv t ) dt ds(p), 
Js Jo 

and another integration by parts proves 



= / / u t vdtds(p)= / / uv t dtds(p). 
Js Jo Js Jo 

Using this and one integration by parts in (JT5j) establishes (flBI) . which completes the proof 
of Theorem [5] for n = 2. The extension to higher (even) dimensions follows almost word for 
word the proof from [51 Section 4.2], where the trace identities in odd dimensions greater 
than three were proved from the three dimensional case. □ 

Proof of Theorem [4] for n > 2. Reversing the chain of reasoning in (I4~8l) proves (|T3l) in 
L 2 sense from (1151) . Similarly, (fT4j) follows from f|T6|) . However, as both sides are continuous 
functions when / is smooth, the formulas hold pointwise as well. □ 



4 Numerical results 



In the previous sections, we have established several exact inversion formulas to recover a 
function / supported in a closed disc D from either its spherical means M. f or the trace W / 
of the solution of the wave equation with initial data (/, 0). However, those formulas require 
continuous data, whereas in practical applications only a discrete data set is available. For 
example, in thermoacoustic tomography (see Figure [T]) only a finite number of positions of 
the line detectors and finite number of samples in time are feasible. In this section we derive 
discrete filtered back-projection (FBP) algorithms with linear interpolation in dimension two 
and present some numerical results. 

The derived FBP algorithms are numerical implementations of discretized versions of 
(SJ-dZj) and f|T2l - f|T4l) and the derivation of any of them follows the same line. We shall focus 
on the implementation of (jSJ), assuming uniformly sampled discrete data 

F Km ;= {M f) (pfc; rm) ; {K m) 6 {0, . . . , N v } X {0, ... , N r } , (49) 
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where p k := Rq, (cos(kh ip ) , sin(kh v )) , r m := mh r , h v := 2n/(N v + 1) and h r := 2Rq/N t . In 
order to motivate the derivation of a discrete FBP algorithm based on (JSj), we introduce the 
differential operator T> := d r rd r and the integral operator 

X : G °°(S x [0, 2Rq)) ^C°°{S x [0, 2Rq)) 

f 2R o _ (50) 

(IG)(p,r) := / G(p, r) log | r 2 — r 2 1 dr 
Jo 

which both act in the second component, and the so called back-projection operator 
B-.C^iSx [0,2R )) -^C°°(D) 

(BG){x) := J G{p,\x-p\)ds(p) ^ 

1 f 2n 

G(p((p),\x-p(<p)\)d(p , 



^JO 

where p(<p) ■= Ro(cos(p,smip). Therefore, we can rewrite (j5J) as 

f = (B!V)(Mf). (52) 

In the numerical implementation the operators B, X, and T> in (|52|) are replaced with finite 
dimensional approximation B, I and D (as described below) and (1521) is approximated by 

j{x i )^f :=(BIDF)% 2G{0,...,iV} 2 . (53) 

Here F := (F k > m \ m with F fc ' m defined by gHD, x i := -(R ,R ) + ih x with z = (i l9 i 2 ) e 
{0,iV} 2 and /i^ := 2Rq/N. In the following S^r and S x denote the sampling operators that 
map G E C°°(S x [0,2^]) and / e G°°(T7) onto its sampks, S^ r G := (G(p fc , r m )) fc>m and 
S x / := f := (f(x' l ))i, where we set f(x l ) := if x % D. Moreover, | • ^ denotes the 
maximum norm on either R,( A V+ 1 ) x ( iV >-+ 1 ) or R,(^+i) x (JV+i) _ 



1. The operator T> can be written as d r + r<9 2 . We approximate d r G with symmetric 
finite differences (G k > m+l - G k ' m ~ 1 ) /(2h r ), d 2 r G by (G k > m+l + G fe ' m_1 - 2G fc ' m ) //i 2 and 
the multiplication operator G \— > rG by point- wise discrete multiplication (G k,m )k, m 
{r m G k,m )k t m- This leads to the discrete approximation 

J) . J^(AV+I)x(^r+1) ^(AT^+^xCiVr+l) 

F ^ (^(DG) fc ' m := ^- ^(m+ i)G fe ' m+1 + (m - ^G^™" 1 - 2mG fe,m j J ^ ^ 

where we set G k, ~ 1 := G k,Nr+1 := 0. The approximation of d r with symmetric finite 
differences is of second order and therefore | (S V)f . V — D S ¥ , )7 .)G| 00 < Gi/i 2 for some 
constant Gi which does not depend on h r . 
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2. Next we define a second order approximation to the integral operator X. This is 
done by replacing G(p k , •) in (15U1) by the piecewise linear spline T k [G] : [0, 2R Q ] — > R 
interpolating G at the nodes r m . More precisely, 

I ' R i ( A V+ 1 )x( Ar >-+ 1 ) _> R^V+-l)x(JV r +l) . q ^ ((lG) k,m ) 

is defined by 



and 



T fc [G](r) := G k ' m + !—ll(G k ' m+1 - G k ' m ) , r e [ r m ,r rn+1 ] (55) 



r2Ro 

(IF) fc ' m := / T fc [G](r)log|r 2 - (r m f\dr 
Jo 

N r -1 / ~ r m '+ 1 

G fc ' m ' / log|r 2 - (r m ) 2 |dr 



(56) 



m'=0 



Wr- 1 f~,k,m'+l _ rik,m' ( pr m ' +1 
+ ^ — — / (r - r m ') log |r 2 - (r m ) 2 |<ir 

For an efficient and accurate numerical implementation it is crucial that the integrals 
in (1551) are evaluated analytically. In fact, by straight forward computation it can be 
verified that 



N r -1 , JVr-1 



(IF 



^ ^ gin £fk,m' _|_ ^ ^ ^^yfe,m'+l Qk,m' 



m'=0 

r m'+l 

a™, : = ( r _ r m ) W |r - r m | + (r + r m ) log |r + r m | - 2r 



(57) 



P _ r m ' n m -L ^ 



^r 2 - (r m ) 2 )log|r 2 - (r m ) 2 | - r 2 



j r =r" 
-ro'+l 



Moreover, using the fact that piecewise linear interpolation is of second order [15] and 
that r i— > log |r — (r m ) 2 | is integrable, it can be readily verified that the approximation 
error satisfies \(S<p ir T — I S^Gl^ < C 2 /i 2 with some constant C 2 independent of h r . 

3. Finally, we define a second order approximation to the back-projection (15T|) . The 
discrete back-projection operator B : R( A V+ 1 ) x ( iV >-+ 1 ) _ > R,(^+i)x(jV+i) ig obtained by 
approximating (15T1) with the trapezoidal rule and piecewise linear interpolation (1551) in 
the second variable, 

(BGr:=— — -^T fe [G](|x l -/|), x'gD, (58) 

and setting (B G) 1 := for x 1 (jL D. It is well known (T5J that both linear interpolation 
in r and the trapezoidal rule in ip are second order approximations and therefore 
| (Sjb i3 — B S (/ p jr )G| 00 < C3 max {/i 2 , ft, 2 } for some constant C3. 
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Algorithm 1 Discrete FBP algorithm with linear interpolation for reconstruction f using 
data F. 



1 

2 
3 
4 
5 
G 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 



h v <- 2ir/(N v + 1) 

ft, r <— 2Ro/N r > initialization 

for m, vnl — 0, . . . , N r do > Pre-compute kernel 

Calculate a™, , b™, according to d56D 
end for 

for k = 0, . . . , iV^ do > Filtering 

for m = 0, . . . , N r do 

pk,m ^_ ( m + \/2)F k ' m+l + (m - l/2)F fe ' m_1 - 2mF fc - m > Equation (I53j) 

end for 

for m — 0, . . . , N r do 

F fc > m 4- ^^ij a™, F fc > m ' + &m' (F k ' m ' +1 - F k ' m ') /h r > Equation ® 

end for 
end for 

for zi, «2 = 0, . . . , N do > FBP with linear interpolation 

i <- (n,«2) 

for fc = 0, . . . , A^, do 

Find me {0,...,N r -l} with r m < \p k - x*\ < r m+1 

T <_ F fc ' m + (r - r m ) (F fc - m+1 - F fc ' m ) //i r > interpolation ||55|) 

f <- f + T/(N V + 1) > discrete back-projection (|5Bj) 

end for 
end for 



The discrete FBP algorithm is given by ( 1531) with D, I, B defined in (15~4"1) . (1571) . ( 1581) and is 
summarized in Algorithm^ Using f(a?) = (S x BTVF) 1 = (S x /)* and f = (BID S^, r F)*, 
the discretization error \f(x l ) — f l \ can be estimated as 

+ \B(S v>r l-IS tp>r )(VF)\ 00 (59) 
+ \BI(S lp , r V-BS tp , r )(F)\ OQ . 

Using the facts that B and I are bounded by some constant independent of h r and that the 
approximation of T>, I, B with D, I, B are of second order, implies that 

IS../-BIDFU <Cmax{h 2 r ,hl} , (60) 

for some constant C independent of h r , h v . This shows that the derived FBP algorithm has 
second order accuracy (for exact data). 



22 





Figure 2: Imaging phantom and data. Left: Imaging phantom / consisting of several 
characteristic functions and one Gaussian kernel. Right: Simulated data F — A4 f. 

In the numerical implementation, the coefficients in fl5^|) . are pre-computed and stored. 
Therefore the numerical effort of the evaluating (|54| is 0{N^N Lp ). Moreover, (|54"1) requires 
0{N r N Lp ) operations and the discrete FBP 0(N 2 N v ), since for all (N + l) 2 reconstruction 
points x l we have to sum over N v + 1 center locations on S. Hence, assuming N ~ N r and 
iV ~ Ntp, Algorithm 1 requires 0(N 3 ) operations and therefore has the same numerical effort 
as the classical FBP algorithm used in x-ray CT [12] . Analogous to the procedure described 
above, discrete FBP algorithms were derived using equation fll]), ([6]) for inverting Ai and 
( TT2|) for inverting W. 

In the following we present numerical results of our FBP algorithms for reconstruction 
the phantom shown in the left picture in Figure [U consisting of a superposition of char- 
acteristic functions and one Gaussian kernel. We calculated the data Ai f via numerical 
integration and the operator W f = d t V f using (TlTT) . Subsequently we added 5% uni- 
formly distributed noise to Ai f and 10% uniformly distributed noise to W /. The results 
for iV = Np — N r — 300 using the algorithms based on (pEj), <^ and (fl2"|) are depicted in 
Figures [31 H] and All implementations show good results although no explicit regularization 
strategy is incorporated in order to the regularize the involved (mildly) ill posed numerical 
differentiation. In particular, §6§ and f|T2|) appear to be most insensitive to noise. However, 
for noisy data, the accuracy of FBP algorithms can be further improved by incorporating a 
regularizing strategy similar to that used in [8]. The derived identities in this article provide 
the mathematical foundation for further development of FBP algorithms for the inversion 
from spherical means and the inversion of the wave equation. 
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Figure 3: Numerical Reconstruction with Algorithm 1. Top: Reconstructions from 
simulated data. Bottom: Reconstructions from simulated data after adding 5% uniformly 
distributed noise. 
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Figure 4: Numerical reconstruction from spherical means with 5% noise added. 

Top: Reconstruction using (j3J). Bottom: Reconstruction using ([6]). 
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Figure 5: Numerical reconstruction using (I12I) from trace W/ of the solution of 
the wave equation with 10% noise added. 
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